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Abstract 

We  report  results  of  a  molecular  dynamics  simulation  study  of  the  high-temperature  pyrolysis  of 
polyethylene  (PE)  and  the  mechanical  and  thermal  properties  of  the  resultant  carbonaceous  char. 
The  microstructure  of  pyrolyzed  PE  samples  was  monitored  during  the  simulation.  Mechanical 
properties  of  the  resultant  char  were  studied  for  char  samples  with  varied  microstructure  (average 
coordination  number  and  ring-size  distributions)  to  establish  structure -property  relationships 
between  mechanical  (thermal)  response  properties  and  microstructure.  We  found  that  this 
relationship  can  be  established  based  upon  the  random  network  theory  of  amorphous  media,  if 
additional  topological  constraints  due  to  rings  are  taken  into  account.  Thermal  conductivity  of 
pyrolytic  char  is  investigated  for  samples  with  different  microstructures  for  a  wide  range  of 
temperatures.  Similar  to  the  mechanical  response,  the  thermal  response  properties  are  correlated 
to  microstructure.  It  is  shown  that  the  behavior  of  the  thermal  conductivity  can  be  well  described 
by  Einstein ‘s  heat  transfer  theory,  if  micro  structure  effecs  are  taken  into  account. 

Introduction 

The  process  of  phase  transformation  under  ultra-high  temperature  treatment  is  of  interest  from 
both  fundamental  science  and  technology  perspectives.  Among  many  others  realizations  of  this 
phenomenon,  pyrolysis  of  polymeric  materials  has  been  extensively  investigated  to  understand 
the  kinetic  and  microstructural  aspects  of  char  formation  [1,2,3].  For  example,  the  chemical 
processes  involved  in  thermal  decomposition  of  polyethylene  have  been  investigated  in  [1,2]. 
Reverse  Monte  Carlo  modeling  of  char  microstructure,  performed  on  industrial  char  samples,  has 
been  reported  in  Ref.  [3].  It  was  found  that  char  contains  buckled  graphene  sheets.  Additionally, 
ultra-high  temperature  pyrolysis  experiments  of  carbon  nanotube-filled  polymer-matrix 
composites  have  recently  been  performed  and  structures  of  the  resultant  char  samples  have  been 
investigated.  These  studies  have  been  performed  for  polymethyl  methacrylate/carbon  nanotube 
[4]  and  polypropylene/carbon  nanotube  [5,6]  nanocomposites.  For  both  materials,  the  integral 
char  yield  was  found  to  increase  as  compared  with  pristine  polymer  systems.  In  general,  it  is 
understood  by  now  that  the  process  of  pyrolysis  of  macromolecular  systems  (pristine  and 
composites)  is  a  complex  phenomenon,  dependent  on  various  aspects  of  initial  polymer- sample 
microstructure  and  pyrolysis  conditions.  The  resultant  microstructure  of  char  samples  has  also 
been  the  focus  of  investigations  in  the  past.  One  important  aspect  of  Ref.  [3]  should  particularly 
be  emphasized.  In  that  work,  the  fact  that  char  can  include  structural  units  of  buckled  graphitic 
sheets  was  particularly  importance  for  microstructure  studies,  as  the  observed  picture  clearly 
differed  from  random  networks  which  had  been  previously  shown  to  provide  an  adequate 
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description  for  a  wide  range  of  amorphous  solids  [7,8].  In  light  of  the  above,  finding  structure- 
property  relationships  for  disordered  systems  with  complex  microstructures  is  of  great 
importance. 

In  the  present  work,  we  study  the  mechanical  and  thermal  properties  of  char,  obtained  by  high- 
temperature  pyrolysis  of  a  model  PE  sample,  with  the  main  emphasis  on  structure-property 
relationships.  The  molten-phase  carbonaceous  char  was  obtained  by  high-temperature  pyrolysis 
using  methodology  previously  described  in  Ref.  [9].  In  brief,  a  pyrolytic  molecular  dynamics 
simulation  approach  was  devised  and  used  to  study  the  process  of  polyethylene  charring.  The 
approach  includes  the  kinetics  of  temperature-mediated  hydrogen  removal  from  the  system  and 
high-temperature  relaxation  of  the  carbonaceous  char.  For  more  details  of  the  approach,  the 
reader  is  referred  to  Ref.  [9].  Moreover,  we  have  recently  performed  more  detailed  simulation 
studies  of  the  pyrolytic  charring  process,  with  particular  emphasis  on  the  kinetics  of 
dehydrogenation  and  the  resultant  amorphous  carbon  structure  evolution  under  different 
conditions.  The  major  focus  of  the  present  study  is  the  basics  of  the  kinetics  of  char  formation 
and  the  mechanical  and  thermal  properties  of  solid-state  char  samples.  To  obtain  amorphous 
char  samples  with  different  microstructures,  the  high-temperature  molten-state  char  sample  was 
quenched  at  different  cooling  rates  and/or  relaxation  times.  A  total  of  five  carbonaceous  char 
samples  (differing  in  micro  structure)  was  obtained  and  investigated  for  mechanical  and  thermal 
response  properties. 


Kinetics  of  pyrolysis 

An  initial  simulation  system  consisted  of  amorphous  PE  cell  with  dimensions  along  the  v-  and  y- 
axis  equal  to  22  A  and  along  the  z-axis  equal  to  87  A.  Periodic  boundary  conditions  were  applied 
to  the  system  in  all  three  directions.  The  system  contained  40  PE  chains,  with  the  total  number  of 
atoms  equal  to  4880,  out  of  which  1600  were  carbon  atoms  and  3280  were  hydrogen  atoms.  The 
Brenner  reactive,  bond-order  potential  (REBO)  was  used  to  describe  interactions  between 
hydrogen  and  carbon  atoms  [10].  In  the  past,  this  potential  has  been  successfully  used  to  model 
dynamic  breaking  and  formation  of  chemical  bonds  under  broad  range  of  conditions. 
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Figure  1.  The  color-coded  map  shows  dynamical  evolution  of  coordination  numbers  of  carbon 
atoms  in  PE  during  pyrolysis  at  T=5000  K.  The  nucleation  of  char  and  the  collapse  of  the 
original  structure  are  illustrated  by  corresponding  changes  in  the  average  coordination  numbers. 
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Figure  2.  The  collapse  of  dehydrogenated  carbon  network  during  cooling:  a)  snapshot  of  the 
initial  simulation  system;  b)  snapshot  of  the  simulation  system  after  t  =90  ps;  c)  the  collapsed 
carbon  network  visualized  at  t  =  200  ps. 

The  dynamics  of  char  formation  is  realized  through  C-H  bond  breaking  and  the  removal  of  free 
hydrogen  atoms  from  the  system.  A  proper  scheme  for  removal  of  hydrogen  atoms  is  essential 
for  formation  of  a  cross-linked  carbon  networks  and  collapse  of  the  dehydrogenated  system  into 
a  dense  (highly  coordinated)  carbonaceous  char  structure.  In  real  polymer  systems,  free  hydrogen 
atoms  will  diffuse  out  of  the  system  at  elevated  temperatures  through  a  free  surface  or  some 
other  interface.  To  accelerate  the  process  of  hydrogen  removal  (which  is  essential  to  simulate 
charring  within  MD  time-scales),  hydrogen  atoms  were  removed  from  the  system  after  they 
became  free  of  the  polymer  during  the  course  of  dynamics.  A  hydrogen  atom  was  considered 
free  if  it  was  outside  the  range  of  the  Brenner  potential.  In  Fig.  1 ,  a  color-coded  map  showing  the 
dynamical  evolution  of  the  average  carbon  coordination  number  is  presented  for  the  PE  sample 
pyrolysed  at  T=5000  K.  The  simulation  cell  was  divided  into  slices  of  thickness  ~1  A  along  the 
z-direction  (parallel  to  the  x-y  plane).  The  average  coordination  numbers  for  C  atoms  were 
computed  in  each  bin  and  plotted  as  a  function  of  time.  In  the  figure,  the  average  coordination  of 
1  is  depicted  in  red  color,  while  an  average  coordination  of  2,  3,  and  4  is  shown  in  green,  purple, 
and  blue,  respectively.  As  can  be  observed  in  the  figure,  at  t  =  0,  the  average  coordination 
numbers  are  either  one  or  two,  thus  being  consistent  with  the  original  PE  microstructure.  At 
larger  simulation  times,  the  average  coordination  numbers  changed  from  a  mixture  of  one  and 
two  (red  and  green)  to  pockets  of  threefold  coordination  around  t  =  50  ps.  The  appearance  of 
threefold-coordination  regions  indicates  the  formation  of  graphitic  cross-links  in  the  structure. 
These  cross-links  act  as  nucleation  centers,  which  grow  and  then  merge  forming  larger,  more 
highly  coordinated  structures.  The  interiors  of  these  regions  have  mixed  threefold  and  fourfold 
coordination  numbers,  indicating  smaller,  diamond-like  structures  in  the  core.  By  t  =  80  ps, 
much  of  the  simulation  cell  has  mixed  threefold  or  fourfold  coordination,  but  the  structure  is  still 
distributed  along  the  entire  cell.  By  /  =  100  ps,  however,  the  inner  core  of  the  threefold 
coordinated  regions  starts  to  nucleate  fourfold  regions  and  the  rest  of  the  structure  starts  to  othe 
final  evolved  structure,  thus,  has  highest  coordination  values  three  or  four  in  the  interior  of  the 
char  with  smaller  values  of  3  and  2  observed  toward  the  surface.  Since  the  system  has  collapsed 
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to  less  than  half  the  original  volume  of  the  simulation  cell,  the  density  of  the  carbon  structure  has 
increased  proportionally  to  an  approximate  value  of  1.67  g/cm3  (initial  density  being  0.89 
g/cm3).  One  can  also  observe  that  dynamical  processes  continue  after  the  system  has  collapsed, 
with  small  carbon  fragments  continuing  to  be  ejected  and  absorbed. 

Char  properties:  Young’s  modulus 

The  mechanical  response  properties  of  carbonaceous  char  samples  were  studied  for  five  samples 
with  different  microstructural  characteristics  (i.e.,  average  coordination  numbers  and  ring 
statistics).  To  investigate  the  mechanical  response  properties,  dynamical  loadings  were 
performed  on  each  sample  by  rescaling  the  coordinates  by  a  factor  of  0.001  along  the  sample’s 
largest  dimension  (z-axis).  The  applied  tensile  (or  compressive)  load  corresponds  to  strain  value 
of  -0.1%.  To  model  the  adiabatic  loading,  the  char  samples  after  each  loading  step  were  first 
relaxed  at  elevated  temperatures  for  no  less  than  5000  MD  steps  and,  subsequently,  were  cooled 
down  to  T=5  K.  The  total  energy  of  the  loaded  system,  Es,  was  computed  after  each  loading  step. 
The  elastic  energy  due  to  tensile  (or  compressive)  loading  was  computed  as  Es-E0,  where  E0  is 
the  total  energy  of  the  unloaded  sample.  The  step-wise  loading  was  continued  until  the  strain 
reached  the  value  of  -15%.  The  strain  energy  was  found  to  be  a  quadratic  function  of  the  applied 
strain  for  all  considered  samples.  The  Young's  modulus  was  obtained  by  fitting  a  quadratic  form 
to  the  simulation  data. 


To  analyze  the  obtained  data,  we  invoked  the  random  network  theory  of  disordered  (glassy) 
solids.  The  idea  of  representing  disordered  solids  as  random  networks  was  introduced  by 
Zachariansen  in  as  early  as  1932,  based  upon  X-ray  diffraction  experiments  performed  on  SiC>4 
glasses.  Subsequently,  this  concept  has  been  developed  into  a  rigorous  mathematical  theory  of 
network  elasticity  (by  Phillips  [7]  and  Thorpe  [8]).  In  the  following,  we  employ  the  framework 
of  theory  of  random  network  continuous  deformation  to  interpret  our  simulation  results.  First,  let 
us  briefly  overview  the  basics  of  the  theory  [7,8].  We  consider  a  random  network  of  bonds 
connecting  nodes,  where  each  node  corresponds  to  an  atom  in  glassy  material.  The  major 
quantity  of  interest  is  the  average  coordination  number,  <r>,  defined  as  <r>  =  Yjnrr[jir,  where 
r  is  the  number  of  bonds  and  nr  is  the  number  of  atoms  having  r  bonds.  Dependent  on  network 
connectivity,  a  random  network  can  be  considered  as  a  polymeric  glass  or  amorphous  solid.  The 
transition  between  these  two  regimes  occurs  at  <r>=  rp=  2.4  [8].  These  classification  derives 
from  percolation  theory  and  are  based  on  considering  percolation  of  floppy  (soft)  and  rigid 
regions,  with  amorphous  solids  corresponding  to  a  network  connectivity  above  the  threshold,  i.e., 
<r>  being  greater  than  2.4.  It  should  be  noted  that  the  above  holds  for  random  networks  without 
rings.  As  we  show  below,  ring  corrections  can  introduce  considerable  modifications  to  the  above 
framework.  To  analyze  our  simulation  data  for  the  microstructure  effects  on  the  elastic  response, 
we  use  the  following  expression  for  the  Young's  modulus  dependence  on  <r>: 


r0~  2.4 

where  Yo  is  the  Young's  modulus  of  the  network  with  connectivity  ro  and  rp  is  the  threshold 
value  for  the  percolation  transition,  defined  by  microstructural  properties  of  the  network  [11].  In 
our  analysis,  r0  corresponds  to  100%  sp3  a-C  network,  i.e.,  r0  =  4,  and  Y0  (=823  GPa)  is  the 
Young's  modulus  of  the  100%  sp3  amorphous  carbon. 
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Figure  3.  Young’s  modulus  of  five  carbonaceous  char  samples  is  shown  as  a  function  of  the 
average  coordination  number  (solid  circles).  The  lines  represent  the  dependence  of  the  Young’s 
modulus  on  <r>,  plotted  for  different  values  of  rp. 

As  can  be  observed  in  the  figure,  the  dependence  significantly  differs  from  results  reported  on 
random  networks.  Indeed,  it  is  clear,  that  a  super-linear  dependence  cannot  provide  an  adequate 
description  of  the  observed  behavior.  Consequently,  while  elastic  responses  of  a  purely  random 
network  are  known  to  demonstrate  a  ~  <r>]'5  dependence  with  average  coordination  number,  our 
results  differ  significantly  from  this  type  of  behavior.  One  of  the  key  parameters  of  the  theory  of 
random  networks  is  the  percolation  threshold  rp,  which  separates  the  polymeric  glass  and 
amorphous  glass  phases.  As  the  simulation  results  demonstrate,  our  char  samples  are  not  pure 
random  networks,  as  they  contain  a  non-zero  density  of  rings  of  sizes  varied  from  3  to  8.  The 
effects  of  rings  have  been  previously  considered  in  Ref.  [8].  It  was  shown  that  the  presence  of 
rings  introduces  a  correction  to  the  percolation  threshold  in  the  form  rp  =  2.4  +  5,  where  5  is 
defined  by  densities  of  the  rings  in  the  system.  Following  the  approach  developed  in  Ref.  [8],  we 
found  that  the  correction  factor  5  in  our  char  samples  varies  from  -0.275  to  0.430.  Consequently, 
the  dependence  of  the  elastic  modulus  on  the  coordination  number  in  our  samples  should  include 
the  variance  in  the  parameter  rp.  The  computed  values  of  corrections  (shown  in  Fig.  3)  are  in 
good  agreement  with  the  values  for  this  parameter,  obtained  by  fits  to  the  simulation  data. 
Consequently,  the  elastic  properties  of  amorphous  carbon  are  defined  by  average  coordination 
number  and  a  ring-size  distribution.  The  observed  ring  structures  introduce  a  significant  change 
in  the  mechanical  response  properties  of  char  as  compared  with  purely  random  network  glasses. 

Char  properties:  Thermal  conductivity 

The  thermal  conductivity  of  char  samples  was  investigated  using  the  direct  method  for  thermal 
conductivity  measurement  in  molecular  dynamics.  This  scheme  is  based  upon  the  Fourier  law  of 
heat  conduction  [12].  The  adopted  scheme  is  as  follows.  The  simulation  cell  was  divided  along 
the  z-axis  into  100  equal-size  units.  Two  slabs,  located  at  Lzl 4  and  -Lz/4  with  respect  to  the 
center  of  the  system,  are  used  as  a  heat  source  and  a  heat  sink,  respectively.  The  velocity  of  each 
atom,  residing  in  these  slabs,  was  rescaled  at  each  MD  step  to  increase/decrease  the  kinetic 
energy  by  a  fixed  value  ±Ae,  such  that  the  total  energy  of  the  system  is  conserved.  The  procedure 
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is  continued  until  the  system  is  equilibrated,  with  the  temperature  gradient  reaching  a  constant 
value.  The  thermal  conductivity  is  computed  using  the  following  expression 


(2) 


where  X  is  the  thermal  conductivity  and  Jz  =  A£/(2LxLyAt)  is  the  heat  flux.  The  temperature 
dependence  of  the  thermal  conductivity  of  our  char  sample  with  largest  coordination  number 
(<r>=  3.10)  is  presented  in  Fig.  4.  As  can  be  observed  in  the  figure,  the  thermal  conductivity 
gradually  increases  with  temperature  in  the  studied  range  of  temperatures,  with  behavior  at 
highly  elevated  temperatures  indicating  near  saturation,  compared  with  the  fast  growth  in  the  low 
and  moderate  temperature  regimes.  This  behavior  can  be  understood  in  terms  of  incoherent 
atomic  vibrations.  As  the  temperature  increases,  the  mean  free  path  of  such  vibrations  decreases 
and,  at  very  high  temperatures,  becomes  comparable  with  an  average  inter-atomic  distance.  The 
saturation  value  of  the  thermal  conductivity  can  be  estimated  as  X  =  1/3  Cvv  a,  where  Cv  is  the 
specific  heat,  v  is  the  velocity  of  sound,  and  a  is  the  average  interatomic  distance.  It  should  be 
noted,  however,  that,  in  the  limit  of  very  high  temperatures,  the  behavior  of  the  thermal 
conductivity  can  be  affected  by  enhanced  phonon-phonon  scattering  and  temperature- activated 
microstructural  changes.  These  effects  can  lead  to  decreases  with  temperature  of  the  thermal 
conductivity  in  high  temperature  regimes.  We  observed  structural  relaxation  in  systems  heated 
above  -600  K. 

To  provide  an  analytical  description  of  the  thermal  conductivity,  we  use  Einstein's  theory  of  heat 
transport  [13].  In  the  past,  this  theory  has  been  successfully  applied  to  interpret  experimental  data 
on  the  thermal  conductivity  of  disordered  solids  [13]  and  amorphous  carbon  thin  films  [14]. 


Figure  4.  The  thermal  conductivity  versus  temperature  is  shown  for  char  sample  with  <r>= 
3.10.  Shown  in  the  inset  are  the  behaviors  of  the  thermal  conductivity  with  atomic  density 
(measured  in  1/A3)  for  different  values  of  Young’s  modulus  (decreasing  from  1  to  5). 
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The  approach  is  based  upon  the  assumption  of  diffusive  heat  transfer  by  incoherent  atomic 
vibrations.  The  heat  diffusion  takes  place  between  Einstein  oscillators  on  a  time  scale  of  half  of 
the  period  of  vibration.  Within  this  approach,  the  minimum  value  of  the  thermal  conductivity  can 
be  estimated  as 
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In  the  above  equation,  index  i  =1,2,3  marks  two  transverse  and  one  longitudinal  sound  modes, 
traveling  with  speeds  v,-,  0i=v,  (hlk^)  (6n 2n)113  (i=  1,2,3)  are  the  cut-off  frequencies  of  these  three 
modes  (expressed  in  Kelvin),  and  n  is  the  number  density  of  carbon  atoms  (expressed  in  1/A3). 
The  analytical  dependence  of  the  thermal  conductivity  as  a  function  of  temperature,  obtained 
using  Eq.  (3),  is  shown  in  Fig.  4  as  a  solid  line.  As  can  be  observed  in  the  figure,  the  theoretical 
approach  provides  an  appropriate  qualitative  description  of  the  behavior  by  capturing  the  major 
trends  of  the  behavior.  This  is  particularly  true  for  the  behavior  at  elevated  temperatures  (for 
temperatures  exceeding  100  K).  In  the  low  temperature  regime,  deviations  from  the  model 
prediction  take  place.  These  deviations  are  due  to  substantial  differences  between  classical  MD 
temperatures  and  quantum  temperatures.  In  other  words,  in  the  low-temperature  regime,  quantum 
corrections  become  significant  and  should  be  taken  into  account  to  achieve  a  better  description  of 
the  thermal  conductivity.  The  analytical  form  Eq.  (3)  can  be  used  to  predict  the  thermal 
conductivity  dependence  for  samples  differing  in  microstructure.  In  the  inset  to  Fig.  4,  we  show 
the  thermal  conductivity  dependence  on  atomic  density  for  samples  with  different  elastic 
properties.  We  used  Eq.  (1)  to  relate  the  microstructure  to  elastic  modulus.  Knowing  the  elastic 
properties  of  samples,  it  is  straightforward  to  obtain  the  corresponding  velocities  of  sound.  Note 
that  the  obtained  ranges  of  thermal  conductivity  variations  are  in  good  agreement  with 
previously  reported  experimental  data  [14]. 


Summary  and  conclusions 


In  summary,  MD  simulations  with  Brenner  bond-order  reactive  potential  have  been  employed  to 
study  mechanical  and  thermal  response  properties  of  carbonaceous  char  samples  with  differing 
microstructure.  Hot  char  samples  were  prepared  using  previously  developed  pyrolytic  MD 
method  as  applied  to  polyethylene.  A  set  of  char  samples  with  varied  microstructure  has  been 
generated  using  different  annealing  regimes  applied  to  the  hot  char.  We  used  the  results  on 
microstructure  to  provide  an  interpretation  of  thermal  and  mechanical  behavior  based  upon 
variations  in  char  microstructure.  First,  we  computed  the  mechanical  responses  of  five  char 
samples  (differing  in  microstructure)  and  interpreted  the  obtained  results  in  terms  of  random 
network  theory  of  amorphous  solids  taking  into  account  corrections  due  to  rings.  The  results 
allowed  us  to  predict  the  mechanical  response  as  a  function  of  microstructural  properties.  Using 
the  direct  method  for  thermal  conductivity  measurement  in  molecular  dynamics,  we  investigated 
the  thermal  conductivity  of  char  samples.  We  found  that  the  thermal  conductivity  can  be  well 
described  by  the  "incoherent  vibrations  model",  with  all  the  essential  features  of  its  behavior 
captured  within  the  framework  this  model.  This  includes  the  temperature  dependence  of  thermal 
conductivity  and  its  dependence  on  the  microstructural  characteristics. 
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